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In this work we numerically demonstrate both significant transient (i.e. non- 
modal) linear amplification and sustained nonlinear turbulence in a kinetic plasma 
system with no unstable eigenmodes. The particular system considered is an electro¬ 
static plasma slab with magnetic shear, kinetic electrons and ions, weak collisions, 
and a density gradient, but with no temperature gradient. In contrast to hydrody¬ 
namic examples of non-modal growth and subcritical turbulence, here there is no 
sheared flow in the equilibrium. Significant transient linear amplification is found 
when the magnetic shear and collisionality are weak. It is also demonstrated that 
nonlinear turbulence can be sustained if initialized at sufficient amplitude. We prove 
these two phenomena are related: when sustained turbulence occurs without unsta¬ 
ble eigenmodes, states that are typical of the turbulence must yield transient linear 
amplification of the gyrokinetic free energy. 


I. INTRODUCTION 

Despite the usefulness of linear eigenmode analysis in many contexts, it has been estab¬ 
lished that such analysis can in practice give misleading predictions about the stability of 
various systems [lj . In particular, turbulence may be sustained in a number of fluid systems 
despite the absence of any unstable eigenmodes. The first such systems known were certain 
sheared flow patterns in neutral fluids, such as plane Poiseuillc or plane Couette flow [23-|3]. 
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Later, sustained turbulence was seen in fluid simulations of several plasma systems with no 
absolute linear instabilities [3HE]. More recently, turbulence in the absence of linear insta¬ 
bility has been seen in kinetic plasma simulations involving sheared equilibrium flow pin]. 
These examples demonstrate that linear eigenmode analysis can give deceiving predictions 
for nonlinear behavior, a conclusion supported by other studies of plasma turbulence in 
which linear instability exists, but the fluctuations in saturated nonlinear turbulence are 
different in character from the linear eigenmodes mm- 

For the cases of turbulence without absolute linear stability in the presence of sheared 
flows, for both neutral fluids and plasmas, the phenomenon of transient linear amplification 
appears to be important [D p snnHEj. Also known as non-modal amplification, transient 
amplification is a period of growth before the exponential decay sets in. A necessary (but 
not sufficient) condition for transient amplification is that the linear operator must not have 
a complete orthogonal set of eigenvectors. Analysis of transient amplification has also been 
applied to plasma systems in which linear instability exists but transient linear behavior 
may be more important [131 dHH22]. 

In the present work, we consider the phenomena of transient amplification and turbulence 
without instability for the case of the weakly-collisional plasma drift wave driven by a density 
gradient [231 l2ij. neglecting toroidal effects. This oscillation has been called the “universal 
instability” as it was at first mistakenly thought to be linearly unstable in any plasma 
with a density gradient. The oscillation is electrostatic in nature, and linear instability 
can indeed be found in slab geometry without magnetic shear or temperature gradients. It 
was later established that when any amount of magnetic shear is included in the model, 
no matter how small, the eigenmodes become linearly damped for k y pi -C 1, where k y is 
the perpendicular wavenumber and pi is the ion gyroradius [25H2T] . However, instability 
can persist at k y pi > 1 when the collisionality is very low [28]. When both magnetic shear 
and modest collisions are included [23], the mode becomes linearly stable at all k y pi. It 
is this latter situation that we will consider here. (As an aside, the small-but-finite-shear 
eigenmodes are different from zero-shear modes not only for the universal instability, but 
also for the slab ion-temperature-gradient (ITG) instability. Specifically, in [29], considering 
the slab mode by setting curvature drift effects Ud to 0, the small-shear limit of (19) is 
different from the zero shear result (12).) 

The fluid version of this drift wave system was investigated in [8], using a model with 
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three fluctuating quantities (density, electrostatic potential, and parallel velocity). In this 
model all linear modes were damped, yet nonlinearly, turbulence could be sustained. It was 
pointed out that an important secondary instability in this system is a drift wave rotated 
by 7 t/ 2 in the plane perpendicular to the magnetic field, since the stabilizing effects of 
magnetic shear are reduced in this orientation. In this rotated drift instability, the role of 
the equilibrium density gradient is replaced in a sense by the perturbed density gradient, 
the latter represented in the nonlinear term. However, since the nonlinearity is conservative, 
the rotated drift wave does not inject energy into the system in the way that can be done by 
a perturbation feeding off the density gradient of the equilibrium. Several authors have also 
investigated transient linear effects in fluid drift wave models mmm- However, these 
studies did not include magnetic shear, which we will find to have a very strong effect on 
transient behavior. 

In the following sections, we establish three results concerning the drift wave system 
in slab geometry with magnetic shear, extending previous work by considering a kinetic 
description and exploring the role of transient linear amplification. For these studies, the 
initial-value gyrokinetic code gs2 [30j will be used for numerical computations. First, in 


section |III| we demonstrate that significant transient linear amplification is possible in this 
system, and that the transiently growing fluctuations resemble the zero-shear eigenmodes in 


the dependence of their instantaneous growth rate on wavenumber. Second, in section IV 
we demonstrate that universal-mode turbulence without linear instability, seen previously 
in the fluid model in [8], is also possible in the gyrokinetic model. Finally, we prove there 
is a relationship between transient linear amplification and nonlinear turbulence in this 
system, using similar arguments to the proof for hydrodynamics in [31] and arriving at 
similar conclusions. Specifically, if the physical parameters are such that turbulence can be 
sustained, typical states of the turbulence must be amplified linearly for some amount of 
time, even if these states eventually decay exponentially. Two more precise statements of 
this result and their proof will be given in section [V] 


II. GYROKINETIC MODEL 


We consider slab geometry with magnetic shear, so the magnetic field in Cartesian co¬ 
ordinates ( x',y',z') is B = [e z i + (x' / L s )e y i]B for some constant B and shear length L s . 
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A field-aligned coordinate system is introduced: x — x', y = y' — x'z'/L s , z = z\ so 
B ■ Vx = B ■ Vy = 0 and B • V z = B. In this geometry, we consider the nonlinear 
electrostatic gyrokinetic equation |32j, dropping toroidal effects: 


dh s dh s 1 I'd ($> R s dh s d ($) R dh s 
dt +V{l dz + B dX s dY s dY s dX s 


_<h f 9 (®)r s fMs d( ^)R s 
T/ Ms dt BL n dY s 


-( C s {h s }) Ra 
( £s ~ 2 ) 7,8 _ ’ 


( 1 ) 


together with the quasineutrality condition 



Here s G {i, e} subscripts denote species, <f>(f, r) is the electrostatic potential, h s (t, R s ,S s , //) 
is the nonadiabatic distribution function, £ s = m s v 2 /(2T s ), /i = v\/(2B), and /Ms is 
the leading-order Maxwellian. The scale length of the equilibrium density n is L n = 
—n/(dn/dx), qt = e = — q e is the proton charge, T s is the species temperature, C s is 
the collision operator, and () r and () R denote gyroaverages at fixed particle position 
r = ( x,y,z ) and guiding-center position R s = (X s , Y s , Z s ). The distinction between z 
and Z s may be neglected. Expanding the fluctuating quantities in spatial Fourier modes as 
$ = T,k x ,ky ®k x ,ky(t, z) exp (ik x x+ik y y) and h s = Y, kx) k y h s,k x ,k y (t , z, S 8 , /i) exp (ik x X s +ik y Y s ), 
we may replace the gyro-averages with Bessel functions: (&k x , ky exp(ik x x + ik y y)) R = 
Jos$k x ,k y exp (ik x X s YikyY s ) and {h s , kx ,k y exp(ik x X s + ik y Y s )) r = J 0s h s ,k x ,k y exp (ik x x + ik y y) 
with J 0s = Jo(k±v±/Q s ) and = eB/m s . The perpendicular wavenumber is given by 


k _l = \jk 2 y + (k x - k y z/L s ) 2 . 


(3) 


The system is uniquely specified by the following dimensionless parameters: L s /L n , rj t , r] e , 
T e /Tj, and collisionality v ee L n /vi , where v ee = \Z2nne 4 In A/[( 47 reo) 2 ?Ue' / 2 T e 3//2 ] is the electron- 
electron collision frequency and Vi = \j2Ti/rrii. 

Notice that the parallel coordinate z enters the system (Jl])-(j3]) only through ([3]). Also 
notice k x enters the system Q-Q only through (J3]) and through the nonlinear term in ([Tj). 
Therefore, a linear mode at nonzero k x and nonzero k y satisfies exactly the same equations 
as the corresponding k x = 0 mode, but translated in z by L s k x /k y . As a consequence, the 
linear growth or damping rates of k y 7 ^ 0 modes are independent of k: x , and the eigenmode 
structures are independent of k x up to this translation in 0 . 
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One important property of the system (JT])-(J2]) is the free energy conservation equation 
[331 33] obtained from the f d 3 R s f d 3 v(T s h s / /m s )(. ..) moment of 0: 


dW 

dt 


S(hi, h e ), 


( 4 ) 


where 


is the free energy, 5f s = h s 
from a Maxwellian, and 


W = ^ J d 3 r J 


d 3 v 


Ml 

Ms 


( 5 ) 


Qs^fMs/T s is the total departure of the distribution function 


S(h h h e ) 



d 3 v 


T s h s 

fMs 


c s {h s } 



d 3 v h s 3 )R S T s dn s 
B dY s n s dx 




( 6 ) 


is the net energy source. (To derive 0-0, the equalities f d 3 R s J d 3 v (...) = 

f d 3 R s J d 3 v(.. .) ~ f d 3 r J d 3 v(.. .) = f d 3 r J d 3 v (■ ■ -) r are used.) The first term on the 
right-hand side of (|6]) is negative-definite (by the //-theorem). The second term has the form 
of the particle and heat fluxes multiplied by the density and temperature gradients, and this 
term increases W if the fluxes have the appropriate sign to relax the gradients. Thus, in 
a turbulent steady state where W is statistically constant, equations indicate that 

energy is injected by the equilibrium gradient(s) and dissipated by collisions. Note that the 
operation f d 3 R s used in the derivation of 0 annihilates the nonlinear term in 0, and so 
0-0 apply to both the linear and nonlinear dynamics. 

The system 0-0 may or may not possess absolute linear instability, depending on 
the physical parameters. For the rest of this paper, we take rji = r\ e = 0 so it is certain 
that the ion-temperature-gradient and electron-temperature-gradient (ETG) instabilities are 
suppressed. Even without temperature gradients, the system may still be unstable to the 
universal instability driven by the density gradient. For reference, the growth rates and 
frequencies of this instability in the absence of magnetic shear are given in Appendix [0 As 
detailed in the appendix, some insight into the instability can be obtained by expanding the 
plasma dispersion function Z for the electrons in the small-argument limit, Z(ui/(\k\\\v e )) ~ 
Za /tt for \u}/(k\\v e )\ 1, and by expanding the ion Z function in the large argument limit, 
Z{uj/{\ k\\\vi)) ~ —\k\\\vi/(jj for \uj/(k\\Vi)\ 1 . One then Ends a frequency with positive 

imaginary part oc iy/ir/L n , indicating the instability is due to resonant electrons in the 
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presence of a density gradient. When an arbitrarily small amount of magnetic shear is 
introduced, the global eigenmodes are no longer recognizable as the unsheared slab modes, 
and indeed are absolutely stable for k y pi <C 1 [25112S] • One might reasonably expect weakly 
sheared systems with stable global eigenmodes to exhibit transient linear amplification, 
with the number of e-foldings of amplification tending to infinity as the magnetic shear 
tends to zero. In the following section, we will show this hypothesis to be true, and thereby 
demonstrate the apparently singular limit of zero shear to be a continuous limit of the 
general shear case. 


III. TRANSIENT LINEAR AMPLIFICATION 

Turning now to linear initial-value simulations, we indeed find that significant transient 
linear amplification is possible for certain parameters. Figure [lja illustrates a linear com¬ 
putation for a single (k x ,k y ) in a case of very weak magnetic shear, L s /L n = 300. Other 
parameters used are k y pi = 1 (where pi = Vj/fh), v ee L n /vi = 0.05, nrii = 3672 rn e , T e = Ti, 
and rji = r] e = 0. The results are independent of k x up to a translation in z, as noted 
previously. Collisions are also included to eliminate absolute instability for k y pi > 0.7, as 
discussed in [25], using the collision operators detailed in Ref. [35] . As transient amplifica¬ 
tion generally depends on the choice of norm, here we compare three reasonable quadratic 
norms: W, f gR| < F| 2 , and the maximum of |<f>| 2 along z. The curve for each norm is scaled 
by a constant so the initial amplitude is 1. Amplification by several orders of magnitude is 
observed in all norms before exponential decay sets in. 

Results for the three norms shown are quite similar, though the curve for the W norm 
appears slightly higher than the other curves due to less decay in the initial period tvi/L n < 
20. Indeed, we find all results described or plotted for the remainder of this paper are quite 
similar for all of these norms. Henceforth we restrict our attention to the W norm since it 
is well motivated by the conservation equation Q. 

For the calculation in figure [lj the initial condition used is the standard “noise” initial 
condition of gs2: g s = h s — q s fMs/T s is set to a Maxwellian velocity distribution 

times a random complex number at each grid point in z. (The real and imaginary parts 
of these complex numbers are each taken from a uniform distribution on [—1,1].) While 
techniques exist to calculate the optimal perturbation that gives maximum amplification 
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(BH2HI22], such techniques are numerically challenging for this kinetic model due to the 
high dimensionality of the problem. The noise initial condition is a convenient shortcut, 
ensuring energy is likely to be given to a growing perturbation if one exists. It has been 
shown in [19j [22] that for studies of transient linear amplification, random initial conditions 
are a reasonable proxy for a directly calculated optimal initial condition. 






FIG. I: (Color online) (a) Linear evolution for L s /L n = 300, rp = r] e = 0, u ee = 0.05 Vi/L n , and 
kyPi = 1, illustrating transient amplification in three norms. The initial conditions are detailed 
in the main text, (b) The amount of transient amplification may be defined relative to t = 0 or 
relative to the minimum that occurs shortly thereafter, (c) The potential structure along a field 
line for this example. 

Figure[l]b shows the same computation, zooming in to the initial transient period. Before 
the amplification begins around tv\/L n = 20, there is an initial decay period due to strongly 
damped modes. An amplification factor may be defined relative to either to this minimum 
amplitude, or relative to the amplitude at t — 0. In both cases, we define the amplification 



















factor to be 1 if there is no such transient amplification. 

The transiently growing structures in the presence of magnetic shear are closely related to 
the zero-shear eigenmodes. To demonstrate this relationship, we consider a single (k x ,k y ), 
and define an effective k± at each time for the transient fluctuation by 

J dz k±(z ) \$k x ,k y (t, z)\ 2 


(k±) (t) = 


I dz k±(z) \$ kx ,k y (t,z)[ 


(7) 


Similarly, we define an effective /cm at each time for the transient structure by 


f dku 

k\\\ 

2 

®k x ,k y (t,k\\) 

_? 

\ 

c 

1 

a, 

_?r 

®k x ,k y {t, fc N ) 

2 


( 8 ) 


where <$>k x j. y is the Fourier transform of $>k x ,k v i n z. Using these effective values of /cj_ and 
/c||, the zero-shear dispersion relation may be evaluated. The resulting growth rate can 
then be compared to the instantaneous growth rate of the finite-shear transient structure, 
7 = (2 W(t))~ 1 dW/dt. This comparison is shown in figure [2j In this figure, we display 25 
independent simulations of transients, varying L s /L n between 20 and 300, and plotting the 
evolution of 7 from the time of peak growth rate to the time at which the energy begins to 
decay. In each case, there is very close agreement between the instantaneous growth rate of 
the transient and the ‘predicted’ growth rate from the zero-shear dispersion relation. Thus, 
we can conclude that the transiently growing structures are closely related to the zero-shear 
eigenmodes. Magnetic shear causes k± and k\\ to vary with time as the structure evolves 
(unlike the zero-shear case in which k± and k\\ are constants.) However, at each instant, the 
finite-shear transient grows at precisely the rate expected from the zero-shear instability at 
the average wavenumbers (k±) and (|fc|||)- For the simulations in figure [ 2 J k y pi = 2 and 
Vee 0.05 V{/ L n . 


For the range of k\\ and k± values computed for the transient structures of figure [2j the 
zero-shear dispersion relation is well approximated by the approximate analytic dispersion 


relation (A2). Thus, we can also conclude that the transient structures, like the zero-shear 
eigenmodes, are driven by the combination of resonant electrons and the density gradient. 

Transient amplification is exponentially reduced by increasing magnetic shear, as illus¬ 
trated in figure [3j For this figure, k y pi is fixed at 2, and other parameters are the same as 
in figure [I] (We choose k y pi = 2 here since wavenumbers near this value will turn out to 
be the most amplified, as we will find shortly.) There is some scatter in figure [3j associated 
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Growth rate of finite-shear transient: (L n /vi)(2W) 1 dW/dt 


FIG. 2: (Color online) Although a change from zero to small finite magnetic shear causes a funda¬ 
mental change in the eigenvalues, (which change from unstable to decaying), the transient finite- 
shear structures remain similar to the zero-shear modes. For as shown here, the instantaneous 
growth rate of the transient finite-shear structure is nearly identical to the growth rate of a zero- 
shear eigenmode at the dominant k± and /c| values 0-0- Each trajectory shown is an independent 
linear initial-value calculation, with L s /L n varying from 20-300 between simulations. The growth 
rate is plotted from the time of maximum growth rate (at the + symbols) until decay begins. 

with the random initial condition, but the exponential trend with L s /L n is clear. Notice 
that amplification only appears possible if L s /L n > 20. Values of L s /L n far in excess of 
20 are obtained in the edge of tokamak experiments. (For a high aspect ratio tokamak, 
L s ~ qR/s where q is the safety factor, R is the major radius, s = (r/q)dq/dr, and r is 
the minor radius.) Figure [3}b demonstrates that the time during which amplification occurs 
varies approximately linearly with L s /L n . 

Figure [4] shows the dependence of the amplification factor on k y at fixed L s /L n = 200. 
Significant amplification can be observed for a range of k y surrounding 1/p*. Results are 
shown for two values of collisionality, v ee L n /vi = 0.1 and 0.5. (When the collisionality is 
much smaller than these values, absolute linear instability is present [IB]-) It can be seen that 
increasing collisionality tends to limit amplification at larger k y while having much less effect 
at smaller k y . This kinetic problem is similar to conventional hydrodynamics in that for both 
cases, non-modal amplification decreases with increasing collisional dissipation (associated 



10 


with viscosity in the Navier-Stokes equation for the hydrodynamic case.) However, the 
kinetic problem here is rather different from the hydrodynamic case in the sensitivity of 
the transient amplification to the geometric parameter L s /L n , which is not associated with 
entropy-producing collisions. 
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FIG. 3: (Color online) The (a) number of e-foldings of transient amplification and (b) duration of 
the amplification both scale with the magnetic shear length scale L s , consistent with the appearance 
of absolute instability in the limit of vanishing shear (L s —> oo). Parameters used are k y pi = 2 and 
Vee — 0.05?y/ L n . 



FIG. 4: (Color online) Transient amplification can be observed for a range of k y surrounding 1/pi. 
For this figure, L s /L n = 200. The two collisionalities shown are v ee L n /vi = 0.1 and 0.5. 
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IV. TURBULENCE 


We now demonstrate that when the nonlinear term in (JTj) is retained, sustained turbulence 
is possible in this system. For the demonstration here we consider the parameters L s /L n = 
20, T e = Ti, rrii = 3672 m e , ry = r} e = 0, and u ee L n /vi = 0.05. We find this point in 
parameter-space is approximately on the edge of where sustained turbulence is possible: if 
the magnetic shear or collisionality is increased much beyond these values, turbulence cannot 
be sustained indefinitely (though it may persist for a finite time.) Note that L s /L n = 20 is 
also approximately the maximum shear in figure [3] for which transient linear amplification is 
observed for random initial conditions. (While absolute linear instability was found in [2H] 
for L s /L n > 17 in the limit of vanishing collisionality, here the finite collisionality is enough 
to suppress this linear instability.) Before presenting nonlinear simulations, we first illustrate 
the linear behavior of the system using all the physical and numerical parameters that will 
be used for the nonlinear runs. Figure [5] shows the time evolution of W for all the dealiased 
k y modes that will be used for nonlinear simulations, including twist-and-shift boundary 
conditions, but with the nonlinearity turned off, and using the random initial conditions 


described in section III, More precisely, the quantities plotted are Wk x ,k v , summed over the 


k x that link to k x = 0, where 

Wk x ,k y — L x Ly 'y ( J dz J 


d A v 


T ,, 


2/m, 


| Sf a 


(9) 


are the 2D Fourier contributions to W — ky W kxt k y - Here, L x and L y are the sizes of 
the simulation domain in x and y, and 5f s = J2k x k y $fs,k x ,k y exp {ik x x + ik y y). Curves for 
the other k x in the simulation look similar to those in figure [5j as expected from the z- 
translation symmetry of the system 0-@; see discussion following (|3j) . The amplitude of 
each of these curves decreases for large t, confirming there is no instability for these physical 
and numerical parameters. The decrease is mostly monotonic, but there is typically a small 
amount of transient growth by a factor < 2x for k y pi between 1-2. Thus, this point in 
parameter space is on the edge of the region in which transient linear amplification occurs, 
as anticipated from figure [3j 

Throughout this section, the numerical parameters used were as follows: 12x2 grid points 
in pitch angle and sgn(u||), 12 grid points in energy, 24 grid points per cell in z (with cells at 
k y = jAky linked if separated in k x by 4jAfc x ), CFL number 0.5, 256 x 64 modes in x and 
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y before dealiasing, and box size (L x , L y , L z ) = (31.4pj, 62.8 pi, 160 L n ). It was verified that 
both the linear and nonlinear results reported in this section were not significantly altered 
under factor-of-2 changes in any of these numerical parameters. 

One noteworthy linear property apparent in figure [5] is that damping of modes with 
k y = 0 is very weak, nearly invisible on the scale of the figure. Linear modes with k y = 0 
and d/dz = 0 are damped only by collisions. (In the absence of any collisions, any h e and hi 
that depend only on £ s and /i are steady-state solutions of the linearized system.) Collisions 
introduce damping of these k y = 0 modes, with a damping rate that increases with k^ due 
to the classical transport terms in the gyroaveraged collision operator. 



FIG. 5: (Color online) Linear dynamics for the parameters used in the nonlinear simulations of 
figure [6j The large-t behavior demonstrates there is no absolute linear instability. A small amount 
of transient amplification can be seen. 

Finally, figure [6] shows two nonlinear simulations for the same parameters as figure [5j 
differing only in the initial amplitude of the k y = 0.1 part of the perturbation. For sufficiently 
large initial amplitude, turbulence is sustained, demonstrating nonlinear instability. 

A variety of strategies can be used to initiate the nonlinear instability in numerical 
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simulations. Large-amplitude noise can be used, but this method is numerically inconvenient 
(at least in gs2) due to the rapid drop in amplitude of strongly damped modes over the first 
few time steps, which causes significant changes in the CFL-stable time step associated with 
the perpendicular ExB speed, thereby requiring recomputation of the implicit time-advance 
operator. Another method for initiating the nonlinear instability is to begin with a sufficient 
temperature gradient to bring about linear instability, and then to reduce the temperature 
gradient once the turbulence reaches a sufficient amplitude 0. A third effective method, the 
one used in figure [6j is to initialize with one of the least-damped linear eigenmodes scaled 
up to appropriate amplitude, with noise in the other modes. We find this method to be 
particularly convenient since variation in the stable time step is minimized. 

As another demonstration that the system displays nonlinear instability without linear 
instability, figure [7] shows the result of turning off the nonlinear term after a period of 
saturated turbulence. The nonlinear term is turned off at time tVi/L n = 2064, shown as 
a vertical dotted line. All modes eventually decay. However some modes display a modest 
amount (< 3x) of transient amplification before decaying. 

The particle and heat fluxes for a long nonlinear simulation at these parameters are 
plotted in figure [8] As the time average of S in (|4])-([6]) must vanish, and since the collision 
term is negative-definite, there must be an average particle flux in the appropriate direction 
to relax the equilibrium density gradient, and indeed we find such a flux. Although there is 
no equilibrium temperature gradient, nonzero average fluxes of both ion and electron heat 
are found. Both heat fluxes have the same sign as the particle flux, meaning heat flows in 
the -Vn direction. 

The fluxes in figure [8] are normalized by gyro-Bohm values defined using the length L n 
and including the 2 in the thermal speed, e.g. the particle flux is normalized by (pj/L n ) 2 my. 
Converting to dimensional units for a deuterium plasma, we obtain the following estimate 
for the turbulent particle diffusivity: D ~ (0.4 m 2 /s) (T/500 eV) 3 ^ 2 (B /2 T) -2 (L n /1 cm) -1 . 
This magnitude of diffusion coefficient may be large enough to be of interest experimentally. 
A larger diffusion coefficient is obtained if L s /L n is increased. 
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Normalized time tVi/L n 


FIG. 6: (Color online) Simulations for exactly the same physical and numerical parameters as 
the decaying linear simulation in figure ([5]), but now including nonlinearity and varying the initial 
amplitude of the k y = 0.1 component, (a) When the initial amplitude is sufficient, turbulence 
grows and is sustained, whereas (b) for small initial amplitude, all modes decay as they do in a 
linear simulation. 

V. TURBULENCE IMPLIES TRANSIENT AMPLIFICATION OF TYPICAL 

STATES 


It seems plausible that the previous two sections are related, that a connection exists 
between transient linear amplification and sustained nonlinear turbulence. Let us now prove 
this relationship rigorously. Here we present two such arguments, both of which yield the 
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FIG. 7: (Color online) Demonstration of the concepts in section [VJ Curves show free energy in 
the 2D modes with k x = 0 (including intervals linked by the twist-and-shift parallel boundary 
condition) for various k y . After a period of saturated turbulence, the nonlinear term is turned off 
at normalized time 2064, denoted by the vertical dashed line. All modes eventually decay, but 
transient amplification is apparent first for some of the 2D modes. Such amplification must occur 
in at least one 2D mode in such a linear “restart” from steady turbulence, as proved in section |V} 

conclusion that non-modal amplification is a necessary condition for sustained turbulence. 
One argument for the case of hydrodynamics was given in EH and section 4 of ra. and 
below we adapt this reasoning to gyrokinetics. Before we do so, we give an alternative 
argument that relics less on linear algebra. For this discussion, it is useful to define a “2D 
mode” to be any structure with given k x and k y , but with no restriction on the z dependence. 

Statement 1: Consider a turbulent solution of the nonlinear gyrokinetic system 0-0. 
i.e. a solution for which W does not decay to zero at large t. Choose any time f 0 at which the 
free energy satisfies dW/dt > 0. Suppose at t 0 we turn off the nonlinear term and re-evolve 
the system in time, as in figure [7J Then at least one 2D mode must grow in the W norm 
for some time following to , i.e. there must be transient linear amplification. 






















16 



FIG. 8: (Color online) Fluxes for the simulation of Figure ([6]). a with v ee L n /vi = 0.05 and L s /L n = 
20, approximately the maximum shear at which turbulence can be sustained at this collisionality. 
Even in the absence of equilibrium temperature gradients, nonzero average ion and electron heat 
fluxes are obtained. 

Proof: Recall that Q is a moment of both the nonlinear and linear dynamical equations. 
When the nonlinear term is turned off at to, S(hi, h e ) does not immediately change, because 
the system state does not immediately change. Therefore, dW/dt must vary continuously 
when the nonlinear term is turned off, and so dW/dt is positive for some period following 
to under the action of the linear dynamics. Recall W is a sum of positive contributions 
(|9j) from each 2D mode: W = ^2 kx k ^ Wk x ,k y - Therefore dWk x ,k y /dt must be positive in the 
linear dynamics for at least one 2D mode, so this mode must be growing transiently in the 
W norm. 

Roughly speaking, the upshot of this argument is that we should not be surprised that 
there is transient amplification after the nonlinear term is turned off in figure ([?]). Such 
amplification must occur. Put another way, a necessary condition for turbulence (in the 
absence of linear instability) is that typical states of the turbulence can be transiently 
amplified. Furthermore, the energy growth in the 2D modes that are amplified must be 
sufficient to balance the energy decay from the damped 2D modes. 

Next, we present a second argument for the necessity of transient amplification. This 
argument is based on linear algebra and closely follows EH- For this second approach, we 
let h denote a column vector that specifies the state of the system { h t , h e } in some convenient 
basis. For example, h could represent the array of values specifying the distribution function 
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in a code. 

Statement 2: Given a set of physical parameters for which the nonlinear gyrokinetic 
system (JT|)-({2]) permits statistically steady turbulence, there exists at least one state h 0 such 
that (1) dW/dt > 0 if h = h 0 , and (2) the squared projection of the turbulent state h(t) 
onto ho, defined below, has a non-zero average. 

Proof: Let L denote the time-independent linear operator associated with 0-([2| in the 
absence of the nonlinear term, so the linearized ([Tj)-(J2]) can be written 

dh/dt = Lh. (10) 


(Here we will take L to be real since the equations ([!])-([2]) are real, though the following 
argument may be generalized to complex L by replacing the transpose operation with Her- 
rnitian conjugation.) Also let E be the time-independent operator that defines the free 
energy norm: 

W = h T Eh, (11) 


where T denotes the transpose. Note that E is real and symmetric ( E = E T ). Differentiating 

dW/dt = h T [(ELf + EL] h, (12) 


which is equivalent to Q. As discussed previously, the presence of the nonlinear term does 
not alter Q (since the nonlinear term is annihilated in the integration over the perpendicular 
directions X s and 17), so (12) holds for both the linear and nonlinear gyrokinetic equations. 
Introducing a time average (...) over statistically steady turbulence, then 


(h T [(EL) T + EL] h) = 0. 


(13) 


Next, observe that the operator ( EL) T + EL in (12)-(13) is real and symmetric, so 
it possesses a complete set of orthogonal eigenvectors with real eigenvalues. Let uj and 
A j denote these eigenvectors (normalized so ujuk = Sj^) and their associated eigenvalues. 
In the hydrodynamics literature, the eigenvectors Uk are sometimes called instantaneous 
optimals; the Uj with maximum A j is the fluctuation that maximizes the instantaneous 
dW/dt. We may decompose the state vector in the Uj basis: h{t) = Yhj a j(f) u j for some 


amplitudes aj(t). Using this decomposition in (13) gives 


(K WI 2 ) = °- 

3 


(14) 
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Neglecting the uninteresting case in which all terms in the sum (14) are zero, then there must 
be at least one positive eigenvalue A j 0 with an associated nonzero average weight (|aj 0 (f)| 2 ). 
Considering the associated u K in (12), then this state Uj 0 must give positive dW/dt, proving 
statement 2. 

Thus, we have showed in two ways that transient linear amplification is a necessary 
condition for turbulence in the associated nonlinear system. We have also shown that 
the amplitude of the transiently growing state(s) must have a significant amplitude in the 
turbulent state, in the following senses. In the first line of argument, not only must there be 
at least a single growing 2D mode, but there must be enough growing 2D modes of sufficient 
amplitude to balance the negative dWk x ,k y /dt of all the decaying 2D modes. In the second 


argument, in (14) there must not only be at least one state j with an instantaneous growth 
rate Xj > 0, but the amplitudes of states with positive and negative A j must be in balance 
so (dW/dt) = 0. In a statistically steady state, there must be a balance between energy 
input from the instantaneously growing linear states and the instantaneously decaying linear 
states. 

Both of the arguments above relied on the fact that the norm W satisfies the same 
equation Q in both the linear and nonlinear dynamics. Thus, the “nonlinear invariant” W 
is a special norm, and the arguments do not generally apply to other norm-like quantities 
such as f d 3 r |$| 2 in which one might choose to measure transient amplification. (In fact for 
any operator with damped eigenmodes, even if transient amplification exists in one norm, 
another norm always exists in which transient amplification is impossible. The sum of 
squared eigenmode coefficients is such a norm [36].) 

While the arguments here relied on the gyrokinetic conservation law Q, the same reason¬ 
ing may be applied to other models and systems in which an energy conservation equation 
holds. For example, eq (6) in Ref. [8] gives the appropriate conservation law for the fluid 
model in that work, and the first line of (6) gives the associated norm to use in place of W. 


VI. DISCUSSION AND CONCLUSIONS 

The “universal mode” system (slab geometry with a density gradient, electrostatic, no 
sheared equilibrium flow, and weak temperature gradients) is known to have a dramatic 
change in linearly stability as the magnetic shear is raised from zero to small finite values 
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[25H28] . As the eigenvalues are so sensitive to magnetic shear, we might expect traditional 
modal analysis to give a misleading picture of typical behavior, as with other systems that 
have sensitive eigenvalues [T] . Indeed, here we have demonstrated that significant transient 
linear amplification is possible in the drift wave system even when all eigenmodes are decay¬ 
ing. As shown in figure [2j the transient structures grow at the zero-shear growth rate, when 
evaluated at the dominant k± and k\\ at that instant. We have demonstrated that the amount 
and duration of amplification scale with the magnetic shear length L s , so these quantities 
naturally become unbounded in the limit of vanishing shear. The non-modal amplification 
also decreases with increasing collisionality, particularly at the shortest perpendicular wave¬ 
lengths, and the amplification tends to be largest at wavelengths near k y pi ~ 1 — 2. In 
contrast to more famous examples of non-modal amplification from hydrodynamics [3], and 
to some recent examples from plasma physics P GDI QHIIB], the amplification in the plasma 
drift wave system here does not depend on sheared equilibrium flow. 

We have also demonstrated using nonlinear gyrokinetic simulations that sustained turbu¬ 
lence is possible in this system even when all linear modes are damped. This phenomenon 
had been seen previously in fluid models [6] SSj - This phenomenon, clearly visible in figure [6j 
is opposite to the “Dirnits shift ” B3 seen at different physical parameters, in which nonlinear 
fluctuations are suppressed despite the presence of linear instability. 

In this system without unstable eigenvalues, we find that sustained nonlinear turbulence 
and transient linear amplification both become possible around L s /L n > 20 (for the col¬ 
lisionality u ee L n /vi = 0.05 and other parameters considered.) As discussed in section [VJ 
the fact that these two phenomena appear together is not a coincidence. Transient linear 
amplification in the free energy norm W is a necessary condition for statistically steady 
turbulence, and the turbulent state must have a significant projection onto states that are 
linearly growing in the W norm. When the nonlinear term is turned off in a simulation of 
turbulence, as in figure [FJ amplification in W norm must occur in at least one of the 2D 
modes. The amplitude of the instantaneously growing states (instantaneous optimals) must 
be significant, enough for the associated input of free energy to balance decay of energy from 
the states that are decaying. 

Considering the nonlinear and linear results together, the existence of transient ampli¬ 
fication and turbulence evidently depends on whether magnetic shear is sufficiently weak 
that un-sheared slab mode dynamics can persist. This geometric property can be controlled 
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independently of the energy source of the turbulence (i.e. the density gradient), in contrast 
to turbulence driven by sheared flows where the nonuniformity of the flow plays two roles 
at once: shaping the global structure of the growing structures and providing energy to the 
turbulence. This feature makes the present system an interesting case study of sub-critical 
turbulence, and we are unaware of other such examples. 

Note that the transient linear process emphasized here and secondary instability are not 
mutually exclusive. Both are required to support turbulence. Some form of secondary 
instability (as discussed in [S]) must be present to transfer energy out of the transiently 
growing linear structures before these structures begin to decay substantially. However, 
nonlinear processes cannot directly contribute to the system’s free energy W, due to the 
annihilation of the nonlinear term in deriving Q. 

Several effects that may be important for realistic plasma turbulence were not included in 
the present work, including toroidicity and trapped particles, and electromagnetic fluctua¬ 
tions. Therefore, further work is required to determine how relevant the nonlinear turbulence 
considered here is for realistic experimental parameters. However, we at least conclude that 
linear stability need not preclude the possibility of significant turbulent particle and energy 
fluxes, and that nonlinear simulations should sometimes be initialized at finite amplitude to 
examine this possibility. 
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Appendix A: Linear stability in the absence of magnetic shear 


For reference, here we present the growth/damping rates and real frequencies of the 
universal instability in the limit of vanishing magnetic shear and vanishing collisions. In the 
absence of shear, we can replace d/dz in ([Tj) with ik\\ for some constant k\\. Using (J2|, and 
after some algebra, the local dispersion relation may then be written 
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Here, ( s = (n/(|fc||| v s ), v s = \j2T s /m s is the thermal speed, Z s = Z(( s ), Z is the plasma 
dispersion function, Y JS = Ij(b s )e~ bs , Ij is a modihed Bessel function, b s = k 2 ± v 2 s / {2kl 2 ), 
and co* = k y T e /(eBL n ). Figure [9] shows the solution of (Al) with largest imaginary part 
for the case T e = T* and the deuterium-electron mass ratio. We take r]i = r] e = 0 so it is 
certain that the ion-temperature-gradient (ITG) and clectron-temperature-gradient (ETG) 
modes are suppressed. Figure [9] was generated by applying nonlinear root-finding to (Al), 
verifying the solution agreed with initial-value gs2 simulations for several values of k\\ and 
k±. It can be seen that absolute linear instability exists for sufficiently small values of |fc||| 
when k±pi < 40. 



|fc|| Vi/u}*\ 
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FIG. 9: (Color online) Solution of the linear dispersion relation in the zero-shear (local) limit, i.e. 


the root of (Al) with largest imaginary part, for rji = r] e = 0, T e = Ti, and rrq = 3672 m e . The (a) 
imaginary and (b) real parts of the frequency are shown. 


Some insight into the density-gradient-driven instability can be obtained by examining 
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the limit |£ e | <C 1 so Z e « iy/n, and Q| S> 1 so Zi « — 1/Q. Considering solutions with 
|a;/a;* | <C 1, one finds 


^*r 0 i 




+ 1 + iy/n 



(A2) 


The imaginary part of u, proportional to Z e ~ iyfir, evidently arises due to the electron 
resonance. 
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